home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / sph_scat.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  110 lines
  1. ; $Id: sph_scat.pro,v 1.3 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5.  
  6. function sph_scat, lon, lat, f, GS = gs, BOUNDS = bounds, NLON=nlon, $
  7.     NLAT = nlat, GOUT = gsout, BOUT = boundsout
  8. ;+
  9. ; NAME:
  10. ;    SPH_SCAT
  11. ;
  12. ; PURPOSE:
  13. ;    Interpolate to a regular grid given scattered samples on the
  14. ;    surface of a sphere.
  15. ; CATEGORY:
  16. ;    Interpolation.
  17. ; CALLING SEQUENCE:
  18. ;    Result = SPH_SCAT(lon, lat, f)
  19. ; INPUTS:
  20. ;    lon = sample longitudes, a vector, in degrees.  lon, lat, and 
  21. ;        f must have the same number of points.
  22. ;    lat = sample latitudes, a vector, in degreees.
  23. ;    f = data values measured at lon and lat.  f(i) = sample value
  24. ;        at lon(i), lat(i).  
  25. ; KEYWORD PARAMETERS:
  26. ;    GS:      If present, GS must be a two-element vector [XS, YS],
  27. ;          where XS is the spacing between grid points in longitude,
  28. ;          and YS is the spacing in latitude. The default is based on
  29. ;          the extents of lon and lat. If the grid starts at longitude
  30. ;          Lonmin and ends at Lonmax, then the default horizontal
  31. ;          spacing is (Lonmax - Lonmin)/(NX-1). YS is computed in the
  32. ;          same way. The default grid size, if neither NX or NY
  33. ;          are specified, is 26 by 26.
  34. ;    BOUNDS:   If present, BOUNDS must be a four element array containing
  35. ;          the grid limits in longitude and latitude of the output grid:
  36. ;          [Lonmin, Latmin, Lonmax, Latmax]. If not specified, the grid
  37. ;          limits are set to the extent of lon and lat.   Warning:
  38. ;          to cover all longitudes, you must directly specify BOUNDS.
  39. ;    NX:       The output grid size in the longitude direction. NX need not
  40. ;            be specified if the size can be inferred from GS and
  41. ;          BOUNDS. The default value is 26.
  42. ;    NY:       The output grid size in the latitude direction. See NX. 
  43. ;    BOUT:      the actual extent of the regular grid, arranged as in
  44. ;          bounds.  An optional output parameter.
  45. ;    GOUT:     The actual grid spacing, a two element optional output array.
  46. ;    
  47. ; OUTPUTS:
  48. ;    Result = regularly interpolated result.
  49. ; COMMON BLOCKS:
  50. ;    None.
  51. ; SIDE EFFECTS:
  52. ;    None.
  53. ; RESTRICTIONS:
  54. ;    Timing. on a Sun SPARCstation LX producing a 36 x 36 output
  55. ;    grid (1296 points), t is ~ .578 + .00368 * N + 2.39e-06 * N^2.
  56. ;    For example:  
  57. ;    N    16    64    256    1024    4096
  58. ;    Time    .7    .8    1.6    6.6    56
  59. ;    Output points are produced at a rate of approximately 2000
  60. ;    points per second.
  61. ;    
  62. ; PROCEDURE:
  63. ;    This routine is a convenience interface to the Spherical gridding
  64. ;    and interpolation provided by TRIANGULATE and TRIGRID.  The
  65. ;    methods are based on the work of Robert Renka, Interpolation of Data
  66. ;    on the Surface of a Sphere, Oak Ridge Natl Lab Technical Paper
  67. ;    CSD-108.  The procedure consists of generating a triangulation of the
  68. ;    scattered data points, estimating the gradients with a local method,
  69. ;    and then constructing a triangle based interpolant of the data and
  70. ;    gradient estimates.  The interpolant is C(1) continuous.
  71. ; EXAMPLE:
  72. ;    Create 50 random longitudes and latitudes, make a function value,
  73. ;    and then interpolate, obtaining a 360 x 360 array of
  74. ;    10 degree by 5 degree resolution that covers the sphere:
  75. ;
  76. ;    lon = randomu(seed, 50) * 360. -180.  ;Make random scattered points
  77. ;    lat = randomu(seed, 50) * 180. -90.
  78. ;    z = sin(lat*!DTOR)        ;Make a function to fit
  79. ;    c = cos(lat*!DTOR)
  80. ;    x = cos(lon*!DTOR) * c
  81. ;    y = sin(lon*!DTOR) * c
  82. ;    f =  sin(x+y) * sin(x*z)    ;The dependent variable
  83. ;  ** Now, given lon, lat, and f, interpolate the data:
  84. ;    result = sph_scat(lon, lat, f, bounds=[0, -90, 350, 85], gs=[10,5])
  85. ;    
  86. ; MODIFICATION HISTORY:
  87. ;    DMS, November, 1994.  Written.
  88. ;-
  89.  
  90. n = n_elements(lon)
  91. if n ne n_elements(lat) or n ne n_elements(f) then $
  92.     message, 'lon, lat, and f must have the same number of elements'
  93. if n le 3 then $
  94.     message, 'Must have at least 3 points'
  95. ;        Construct bounds if necessary
  96. if n_elements(bounds) ne 4 then $
  97.     boundsout = [ min(lon, max=lonmax), min(lat, max=latmax), lonmax, latmax] $
  98. else boundsout = bounds
  99. ;        Get gs, nx, and ny.
  100. if n_elements(gs) ne 2 then begin
  101.     if n_elements(nx) le 0 then nx = 26
  102.     if n_elements(ny) le 0 then ny = 26
  103.     gsout = [boundsout[2]-boundsout[0], boundsout[3]-boundsout[1]] / $
  104.         float([nx-1, ny-1])
  105. endif else gsout = gs
  106. fcopy = f        ;will be rearranged.
  107. TRIANGULATE, 1.0*lon, 1.0*lat, SPHERE=s, tr, FVALUE=fcopy, /DEGREES
  108. return, TRIGRID(fcopy, SPHERE=s, gsout, boundsout, /DEGREES)
  109. end
  110.